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Removing Camera Shake via 
Weighted Fourier Burst Accumulation 

Mauricio Delbracio and Guillermo Sapiro 


Abstract —Numerous recent approaches attempt to remove 
image blur due to camera shake, either with one or multiple 
input images, by explicitly solving an inverse and inherently ill- 
posed deconvolution problem. If the photographer takes a burst 
of images, a modality available in virtually all modern digital 
cameras, we show that it is possible to combine them to get a 
clean sharp version. This is done without explicitly solving any 
blur estimation and subsequent inverse problem. The proposed 
algorithm is strikingly simple: it performs a weighted average 
in the Fourier domain, with weights depending on the Fourier 
spectrum magnitude. The method can be seen as a generalization 
of the align and average procedure, with a weighted average, 
motivated by hand-shake physiology and theoretically supported, 
taking place in the Fourier domain. The method’s rationale 
is that camera shake has a random nature and therefore each 
image in the burst is generally blurred differently. Experiments 
with real camera data, and extensive comparisons, show that 
the proposed Fourier Burst Accumulation (FBA) algorithm 
achieves state-of-the-art results an order of magnitude faster, 
with simplicity for on-board implementation on camera phones. 
Finally, we also present experiments in real high dynamic range 
(HDR) scenes, showing how the method can be straightforwardly 
extended to HDR photography. 

Index Terms —multi-image deblurring, burst fusion, camera 
shake, low light photography, high dynamic range. 


1. Introduction 

O NE of the most challenging experiences in photography 
is taking images in low-light environments. The basic 
principle of photography is the accumulation of photons in 
the sensor during a given exposure time. In general, the more 
photons reach the surface the better the quality of the final 
image, as the photonic noise is reduced. However, this basic 
principle requires the photographed scene to be static and that 
there is no relative motion between the camera and the scene. 
Otherwise, the photons will be accumulated in neighboring 
pixels, generating a loss of sharpness (blur). This problem is 
significant when shooting with hand-held cameras, the most 
popular photography device today, in dim light conditions. 

Under reasonable hypotheses, the camera shake can be 
modeled mathematically as a convolution, 

V = U'k k ( 1 ) 

where v is the noisy blurred observation, u is the latent sharp 
image, k is an unknown blurring kernel and n is additive white 
noise. For this model to be accurate, the camera movement has 
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to be essentially a rotation in its optical axis with negligible 
in-plane rotation, e.g., m. The kernel k results from several 
blur sources: light diffraction due to the finite aperture, out-of- 
focus, light integration in the photo-sensor, and relative motion 
between the camera and the scene during the exposure. To get 
enough photons per pixel in a typical low light scene, the 
camera needs to capture light for a period of tens to hundreds 
of milliseconds. In such a situation (and assuming that the 
scene is static and the user/camera has correctly set the focus), 
the dominant contribution to the blur kernel is the camera 
shake —mostly due to hand tremors. 

Current cameras can take a burst of images, this being 
popular also in camera phones. This has been exploited in 
several approaches for accumulating photons in the different 
images and then forming an image with less noise (mimicking 
a longer exposure time a posteriori, see e.g., IS). However, 
this principle is disturbed if the images in the burst are 
blurred. The classical mathematical formulation as a multi¬ 
image deconvolution, seeks to solve an inverse problem where 
the unknowns are the multiple blurring operators and the 
underlying sharp image. This procedure, although it produces 
good results m, is computationally very expensive, and very 
sensitive to a good estimation of the blurring kernels, an 
extremely challenging task by itself. Furthermore, since the 
inverse problem is ill-posed it relies on priors either, or both, 
for the calculus of the blurs and the latent sharp image. 

Camera shake originated from hand tremor vibrations is 
essentially random O, ID This implies that the move¬ 
ment of the camera in an individual image of the burst is 
independent of the movement in another one. Thus, the blur 
in one frame will be different from the one in another image 
of the burst. 

Our work is built on this basic principle. We present an 
algorithm that aggregates a burst of images (or more than one 
burst for high dynamic range), taking what is less blurred of 
each frame to build an image that is sharper and less noisy than 
all the images in the burst. The algorithm is straightforward 
to implement and conceptually simple. It takes as input a 
series of registered images and computes a weighted average 
of the Fourier coefficients of the images in the burst. Similar 
ideas have been explored by Carrel et al. Q in the context of 
astronomical images, where a sharp clean image is produced 
from a video affected by atmospheric turbulence blur. 

With the availability of accurate gyroscope and accelerom¬ 
eters in, for example, phone cameras, the registration can be 
obtained “for free,” rendering the whole algorithm very effi¬ 
cient for on-board implementation. Indeed, one could envision 
a mode transparent to the user, where every time he/she takes 
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a picture, it is actually a burst or multiple bursts with different 
parameters each. The set is then processed on the fly and 
only the result is saved. Related modes are currently available 
in “permanent open” cameras. The explicit computation of 
the blurring kernel, as commonly done in the literature, is 
completely avoided. This is not only an unimportant hidden 
variable for the task at hand, but as mentioned above, still 
leaves the ill-posed and computationally very expensive task 
of solving the inverse problem. 

Evaluation through synthetic and real experiments shows 
that the flnal image quality is signiflcantly improved. This 
is done without explicitly performing deconvolution, which 
generally introduces artifacts and also a signiflcant overhead. 
Comparison to state-of-the-art multi-image deconvolution al¬ 
gorithms shows that our approach produces similar or better 
results while being orders of magnitude faster and simpler. The 
proposed algorithm does not assume any prior on the latent 
image; exclusively relying on the randomness of hand tremor. 

A preliminary short version of this work was submitted to 
a conference M- The present version incorporates a more 
detailed analysis of the burst aggregation algorithm and its 
implementation. We also introduce a detailed comparison to 
lucky image selection techniques El, where from a series 
of short exposure images, the sharpest ones are selected, 
aligned and averaged into a single frame. Additionally, we 
present experiments in real high dynamic range (hdr) scenes 
showing how the method can be extended to hand-held HDR 
photography. 

The remaining of the paper is organized as follows. Sec¬ 
tion ini discusses the related work and the substantial differ¬ 
ences to what we propose. Section [nl| explains how a burst can 
be combined in the Fourier domain to recover a sharper image, 
while in Section [Iv] we present an empirical analysis of the 
algorithm performance through the simulation of camera shake 
kernels. Section |V] details the algorithm implementation and 
in Section [Vl| we present results of the proposed aggregation 
procedure in real data, comparing the algorithm to state-of-the- 
art multi-image deconvolution methods. Section VI-D presents 
experiments in real high dynamic range (HDR) scenes, show¬ 
ing how the method can be straightforwardly extended to 
HDR photography. Conclusions are Anally summarized in 
Section [ 


II. Related Work 

Removing camera shake blur is one of the most challenging 
problems in image processing. Although in the last decade 
several image restoration algorithms have emerged giving 
outstanding performance, their success is still very dependent 
on the scene. Most image deblurring algorithms cast the 
problem as a deconvolution with either a known (non-blind) 
or an unknown blurring kernel (blind). See e.g., the review 
by Kundur and Hatzinakos cni, where a discussion of the 
most classical methods is presented. 

Single image blind deconvolution. Most blind deconvolution 
algorithms try to estimate the latent image without any other 
input than the noisy blurred image itself. A representative 


work is the one by Fergus et al CD This variational method 
sparked many competitors seeking to combine natural image 
priors, assumptions on the blurring operator, and complex 
optimization frameworks, to simultaneously estimate both the 
blurring kernel and the sharp image (121, (H, (HI, El, lUSl- 
Fergus et al. CD approximated the heavy-tailed distribution 
of the gradient of natural images using a Gaussian mixture. 
In (121, the authors exploited the use of sparse priors for 
both the sharp image and the blurring kernel. Cai et al. ca 
proposed a joint optimization framework, that simultaneously 
maximizes the sparsity of the blur kernel and the sharp image 
in a curvelet and a framelet systems respectively. Krishnan et 
al. (141 introduced as a prior the ratio between the and the 
^2 norms on the high frequencies of an image. This normalized 
sparsity measure gives low cost for the sharp image. In da 
the authors discussed unnatural sparse representations of the 
image that mainly retain edge information. This representation 
is used to estimate the motion kernel. Michael! and Irani da 
recently proposed to use as an image prior the recurrence of 
small natural image patches across different scales. The idea 
is that the cross-scale patch occurrence should be maximal for 
sharp images. 

Several attempt to first estimate the degradation operator 
and then applying a non-blind deconvolution algorithm. For 
instance, CD accelerates the kernel estimation step by using 
fast image Alters for explicitly detecting and restoring strong 
edges in the latent sharp image. Since the blurring kernel has 
typically a very small support, the kernel estimation problem 
is better conditioned than estimating the kernel and the sharp 
image together. In (TSl, (191 , the authors concluded that it is 
better to first solve a maximum a posteriori estimation of the 
kernel than of the latent image and the kernel simultaneously. 
However, even in non-blind deblurring, i.e., when the blurring 
kernels are known, the problem is generally ill-posed, because 
the blur introduces zeros in the frequency domain. Thereby 
avoiding explicit inversion, as here proposed, becomes critical. 

Multi-image blind deconvolution. Two or more input images 
can improve the estimation of both the underlying image 
and the blurring kernels. Rav-Acha and Peleg (20l claimed 
that ''Two motion-blurred images are better than one,” if the 
direction of the blurs are different. In (IH the authors proposed 
to capture two photographs: one having a short exposure time, 
noisy but sharp, and one with a long exposure, blurred but with 
low noise. The two acquisitions are complementary, and the 
sharp one is used to estimate the motion kernel of the blurred 
one. 

Close to our work are papers on multi-image blind de- 
convolution O, (22l, (23l, (24l, (25\ . In (22l the authors 
showed that given multiple observations, the sparsity of the 
image under a tight frame is a good measurement of the 
clearness of the recovered image. Having multiple input im¬ 
ages improves the accuracy of identifying the motion blur 
kernels, reducing the illposedness of the problem. Most of 
these multi-image algorithms introduce cross-blur penalty 
functions between image pairs. However this has the problem 
of growing combinatorially with the number of images in the 
burst. This idea is extended in O using a Bayesian framework 
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for coupling all the unknown blurring kernels and the latent 
image in a unique prior. Although this prior has numerous 
good mathematical properties, its optimization is very slow. 
The algorithm produces very good results but it may take 
several minutes or even hours for a typical burst of 8-10 
images of several megapixels. The very recent work by Park 
and Levoy |[26l relies on an attached gyroscope, now present 
in many phones and tablets, to align all the input images and to 
get an estimation of the blurring kernels. Then, a multi-image 
non-blind deconvolution algorithm is applied. 

By taking a burst of images, the multi-image deconvolution 
problem becomes less ill-posed allowing the use of simpler 
priors. This is explored in l(27ll where the authors adopted a 
total variation prior on the underlying sharp image. 

All these papers propose kernel estimation and to solve 
an inverse problem of image deconvolution. The main 
inconvenience of tackling this problem as a deconvolution, 
on top of the computational burden, is that if the convolution 
model is not accurate or the kernel is not accurately 
estimated, the restored image will contain strong artifacts 
(such as ringing). 

Lucky imaging. A popular technique in astronomical photog¬ 
raphy, known as lucky imaging or lucky exposures, is to take a 
series of thousands of short-exposure images and then select 
and fuse only the sharper ones 1281. Fried 121 showed that 
the probability of getting a sharp lucky short-exposure image 
through turbulence follows a negative exponential. Thus, when 
the captured series or video is sufficiently long, there will exist 
such a frame with high probability. 

Classical selection techniques are based on the brightness 
of the brightest speckle ||28]| . The number of selected frames is 
chosen to optimize the tradeoff between sharpness and signal- 
to-noise ratio required in the application. 

Others propose to measure the local sharpness from the 
norm of the gradient or the image Laplacian (23, if30l . ED, 
1^ . Joshi and Cohen ED engineered a weighting scheme 
to balance noise reduction and sharpness preservation. The 
sharpness is measured through the intensity of the image 
Laplacian. They also proposed a local selectivity weight to 
refiect the fact that more averaging should be done in smooth 
regions. Haro and colleagues (3^ explored similar ideas to 
fusion different acquisitions of painting images. The weights 
for combining the input images rely on a local sharpness 
measure based on the energy of the image gradient. The main 
disadvantage of these approaches is that they only rely on 
sharpness measures (which by the way is not necessarily trivial 
to estimate) and do not profit the fact that camera shake blur 
can be in different directions in different frames. 

Garrel et al. m introduced a selection scheme for astro¬ 
nomic images, based on the relative strength of signal for 
each spatial frequency in the Fourier domain. From a series 
of realistic image simulations, the authors showed that this 
procedure produces images of higher resolution and better 
signal to noise ratio than traditional lucky image fusion 
schemes. This procedure makes a much more efficient use 
of the information contained in each frame. 

Our paper is based on similar ideas but in a different 



Fig. 1. When the camera is set to a burst mode, several photographs are 
captured sequentially. Due to the random nature of hand tremor, the camera 
shake blur is mostly independent from one frame to the other. An image 
consisting of white dots was photographed with a dslr handheld camera 
to depict the camera motion kernels. The kernels are mainly unidimensional 
regular trajectories that are not completely random (perfect random walk) nor 
uniform. 


scenario. The idea is to fuse all the images in the burst with¬ 
out explicitly estimating the blurring kernels and subsequent 
inverse problem approach, but taking the information that is 
less degraded from each image in the burst. The estimation 
of the “less degraded” information is done in a trivial fashion 
as explained next. The entire algorithm is based on physical 
properties of the camera (hand) shake and not on priors or 
assumptions on the image and/or kernel. 

III. Fourier Burst Accumulation 
A. Rationale 

Camera shake originated from hand tremor vibrations has 
undoubtedly a random nature (H, HI, El- The independent 
movement of the photographer hand causes the camera to be 
pushed randomly and unpredictably, generating blurriness in 
the captured image. Figure shows several photographs taken 
with a digital single-lens reflex (DSLR) handheld camera. The 
photographed scene consists of a laptop displaying a black 
image with white dots. The captured picture of the white dots 
illustrates the trace of the camera movement in the image 
plane. If the dots are very small —mimicking Dirac masses— 
their photographs represent the blurring kernels themselves. 
As one can see, the kernels mostly consist of unidimensional 
regular random trajectories. This stochastic behavior will be 
the key ingredient in our proposed approach. 

Let F denote the Fourier Transform and k the Fourier 
Transform of the kernel k. Images are defined in a regular 
grid indexed by the 2D position x and the Fourier domain is 
indexed by the 2D frequency C- Let us assume, without loss of 
generality, that the kernel k due to camera shake is normalized 
such that f k(x)dx = 1. The blurring kernel is nonnegative 
since the integration of incoherent light is always nonnegative. 
This implies that the motion blur does not amplify the Fourier 
spectrum: 

Claim 1. Let k(x) > 0 and f k(x) = 1 . Then, |^(C)| < I 5 VC- 
(Blurring kernels do not amplify the spectrum.) 

Proof. 

Ho 

□ 


J k(x)e^^'^dx < J \k{x)\dx = j k{x)d: 


lx = 1. 
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Fig. 2. Real camera shake kernels were computed using a sharp snapshot captured with a tripod as a reference. The first row shows a crop of each input 
image (frames 1 to 14) and the proposed Fourier Burst Accumulation (from no additional sharpening). As noted before, the kernels are generally regular 
unidimensional trajectories (second row). The four last columns in the second row show the resultant point spread functions (PSF) after the Fourier weighted 
average for different values of p. The kernel due to the Fourier average with p > 0 is closer to a Delta function, showing the success of the method. The two 
bottom rows show the Fourier real and imaginary parts of each blurring kernel (the red box indicates the 7r/2 frequency). The real part is mostly positive 
and significantly larger than the imaginary part, implying that the blurring kernels do not introduce significant phase distortions. This might not be the case 
for large motion kernels, uncommon in standard hand shakes. 


Most modern digital cameras have a burst mode where the 
photographer is allowed to sequently take a series of images, 
one right after the other. Let us assume that the photographer 
takes a sequence of M images of the same scene u, 

V- = ui^ ki rii^ for i = 1,..., M. (2) 

The movement of the camera during any two images of 
the burst will be essentially independent. Thus, the blurring 
kernels ki will be mostly different for different images in the 
burst. Hence, each Fourier frequency of u will be differently 
affected on each frame of the burst. The idea is to reconstruct 
an image whose Fourier spectrum takes for each frequency the 
value having the largest Fourier magnitude in the burst. Since 
a blurring kernel does not amplify the Fourier spectrum (Claim 
1), the reconstructed image picks what is less attenuated, in 
Fourier domain, from each image of the burst. Choosing the 
least attenuated frequencies does not necessarily guarantee 
that those frequencies are the least affected by the blurring 
kernel, as the kernel may introduce changes in the Fourier 
image phase. However, for small motion kernels, the phase 
distortion introduced is small. This is illustrated in Figure 
where several real motion kernels and their Fourier spectra are 
shown. 


B. Fourier Magnitude Weights 

Let p be a non-negative integer, we will call Fourier Burst 
Accumulation (FBAj to the Fourier weighted averaged image. 


M 


,(x) = r 


-1 


Wi{0 = 




(x), 


(3) 


vi=l 


2^7 = 1 




where Vi is the Fourier Transform of the individual burst 
image Vi. The weight Wi := Wi{Q controls the contribution 
of the frequency ( of image Vi to the final reconstruction Up. 
Given for p > 0, the larger the value of |fi 7 (C )|5 the more 
Vi{C) contributes to the average, refiecting what we discussed 
above that the strongest frequency values represent the least 
attenuated u components. 


The integer p controls the aggregation of the images in 
the Fourier domain. If p = 0, the restored image is just the 
arithmetic average of the burst (as standard for example in 
the case of noise only), while if p ^ oo, each reconstructed 
frequency takes the maximum value of that frequency along 
the burst. This is stated in the following claim; the proof is 
straightforward and it is therefore omitted. 

Claim 2. Mean/Max aggregation. Suppose that Vi{C) far i = 
1 ,. .. ,M are such that |^)^l(C)l = \vi 2 {C)\ = ■■■ = I’CiAOl > 
|vi,+i(C)l > l^i^+AOl > ••• > \vij^(C)\andwi{C) is givenby 
0. If p = 0, then Wi{Q = (arithmetic mean pooling), 

while if p ^ oo, then Wi{Q = ^ for i = ii^...^iq and 
Wi(() = 0 otherwise (maximum pooling). 

The Fourier weights only depend on the Fourier magnitude 
and hence they are not sensitive to potential image misalign¬ 
ment. However, when doing the average in the images 
Vi have to be correctly aligned to mitigate Fourier phase 
intermingling and get a sharp aggregation. The images in our 
experiments are aligned using SIFT correspondences and then 
finding the dominant homography between each image in the 
burst and the first one (implementation details are given in 
Section |V]). This pre-alignment step can be done exploiting 
the camera gyroscope and accelerometer data. 

Dealing with noise: The images in the burst are blurry but 
also contaminated with noise. In the ideal case where the input 
images are not contaminated with noise, 0 is reduced to 



as long as \u\ > 0. This is what we would like to have: a 
procedure for weighting more the frequencies that are less at¬ 
tenuated by the different camera shake kernels. Since camera 
shake kernels have typically a small support, of maximum only 
a few tenths of pixels, their Fourier spectra magnitude vary 
smoothly. Thus, \vi\ can be smoothed out before computing 
the weights. This helps to remove noise and also to stabilize 
the weights (see Section [V|). 
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C. Equivalent Point Spread Function 

The aggregation procedure can be seen as the convolution 
of the underlying sharp image u with an average kernel /cfba 
given by the Fourier weights in 0, 

Up=U'k /cfba + n, (5) 

where 

M \ 

■ fcj(C) I (x), (6) 

i=l / 

and n is the weighted average of the input noise. 

The FBA kernel can be seen as the final point spread 
function (PSF) obtained by the aggregation procedure (assum¬ 
ing a perfect alignment). The closer the FBA kernel is to 
a Dirac function, the better the Fourier aggregation works. 
By construction, the average kernel is made from the least 
attenuated frequencies in the burst —given by the Fourier 
weights. However, since arbitrary convolution kernels may 
also introduce phase distortion, there is no guarantee in general 
that this aggregation procedure will lead to an equivalent PSF 
that is closer to a Dirac mass. 

In Figure we show a series of input images and the 
respective motion kernels. The motion kernels were estimated 
using a sharp snapshot, captured with a tripod. Using the sharp 
reference we solve for a blurring kernel ki minimizing the 
least squares distance to the blurred acquisition Vi, namely, 
Wu^Qfi^ki —Vi\ \ (see e.g., 13^ for a similar setup). The figure 
also shows the real and imaginary parts of the Fourier kernels 
spectra. As one can see for most of the kernels the real part 
is mostly positive and considerably larger than the imaginary 
part. This implies that the less attenuated frequencies will not 
present significant phase distortion. This assumption may not 
be accurate for large motion kernels, an unexpected scenario 
in ordinary camera shake. In this example, as p increases 
the equivalent point spread function gets closer to a Dirac 
function. In particular, the FBA kernel for p > 0 attenuates 
significantly less the high frequencies than the regular arith¬ 
metic average (p = 0), leading to a sharper aggregation. 

IV. Fourier Burst Accumulation Analysis 
A. Anatomy of the Fourier Accumulation 

The value of p sets a tradeoff between sharpness and noise 
reduction. Although one would always prefer to get a sharp 
image, due to noise and the unknown Fourier phase shift 
introduced by the aggregation procedure, the resultant image 
would not necessary be better as p —> oo. Figure shows an 
example of the proposed FBA for a burst of 7 images, and the 
amount of contribution of each frame to the final aggregation. 
As p grows, the weights are concentrated in fewer images 
(Figure c) and d)). Also, the weights maps clearly show 
that different Fourier frequencies are recovered from different 
frames (Figure a)). In this example, the high frequency 
content is uniformly taken from all the frames in the burst. 
This produces a strong noise reduction behavior, in addition 
to the sharpening effect. 



(a) Frames crop 1-7 and the Fourier weights wi for p = 11. 



p = 0 p = 3 p = 7 p = 11 p = 20 p = 50 

(b) Fourier Aggregation results for different p values. 



Fig. 3. Weights distribution of the Fourier Burst Aggregation when changing 
the value of p. As p increases, the weights are concentrated in fewer images 
and the aggregated image becomes sharper but also noisier. The difference 
between (c) and (d) lies in the fact that a large weight may not necessarily 
reflect a large contribution to the flnal image, although this is generally the 
case. 

To make explicit the contribution of each frame to the final 
image, the FBA weighted average 0 can be decomposed into 
its contributions: 

M M 

a(C))(x) := 

i=l i = l 

Each term vi is the result of applying the corresponding 
Fourier weights Wi (a filter) to the respective input frame Vi. 
Figure shows each of these terms for a crop of a real burst. 
As the reader may notice, each frame contributes differently. 
None of the frames capture all the structure present in the final 
aggregated image. 

B. Statistical Performance Analysis 

To show the statistical performance of the Fourier weighted 
accumulation, we carried out an empirical analysis applying 
the proposed aggregation with different values of p. We simu¬ 
lated motion kernels following 0, where the (expected value) 
amount of blur is controlled by a parameter related to the 
exposure time Texp. The simulated motion kernels were applied 
to a sharp clean image (ground truth). We also controlled 
the number of images in the burst M and the noise level in 
each frame s. The kernels were generated by simulating the 
random shake of the camera from the power spectral density of 
measured physiological hand tremor mi. All the images were 
aligned by pre-centering the motion kernels before blurring 
the underlying sharp image. Figure shows several different 
realizations for different exposure values Texp. Actually, the 
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Fig. 4. Anatomy of the Fourier aggregation. The first row shows a crop of each input image vi (frames 1 to 9) and the proposed Fourier Burst Accumulation 
for p = 11 (from no additional sharpening). The second row shows the contribution of each frame to the final aggregation vi (rescaled for better 
visualization). The reader can check that the FBA results from the aggregation of different components presented in different frames. This is also confirmed 
by the Fourier weights distribution shown in the bottom row. The bar plot shown on the right indicates the total contribution of each frame to the FBA image. 
None of the input images contain all the structure presented in the final aggregation. 


amount of blur not only depends on the exposure time but also 
on the focal distance, user expertise, and camera dimensions 
and mass Ea. However, for simplicity, all these variables were 
controlled by the single parameter Texp. 

We randomly sampled thousands of different motion ker¬ 
nels and Gaussian noise realizations, and then applied the 
Fourier aggregation procedure for each different configuration 
(Texp,M, s) several times. The empirical mean square error 
(mse) between the ground truth reference and each FBA 
restoration was computed and averaged over hundreds of 
independent realizations. We decomposed the mean square 
error into the bias and variance terms, namely MSE{up) = 
bias(7Xp)^ -b var(7Xp), to help visualize the behavior of the 
algorithm. 

Figure shows the average algorithm performance when 
changing the acquisition conditions. In general, the larger 
the value of p the smaller the bias and the larger the vari¬ 
ance. There is a minimum of the mean square error for 
p G [7,30]. This is reasonable since there exists a tradeoff 
between variance reduction and bias. Although both the bias 
and the variance are affected by the noise level, the qualitative 
performance of the algorithm remains the same. The bias is not 
altered by the number of frames in the burst but the variance 
is reduced as more images are used, implying a gain in the 
expected MSE as more images are used. On the other hand, 
the exposure time mostly affects the bias, being much more 
significant for larger exposures as expected. 

C. Impact of Burst Misalignment 

Misalignment of the burst will certainly have an impact on 
the quality of the aggregated image. For the general case where 
images are noisy but also degraded by anisotropic blur the 
problem of defining a correct alignment is not well defined. 

In what follows we consider that the burst is correctly 
aligned if each Vi satisfies Vi = u 'k ki ^ rii, with the 
blurring kernel ki having vanishing first moment. That is, 
J /c^(x)x(ix = 0. This constraint on the kernel implies that 
the kernel does not drift the image u, so each Vi is aligned 
to u (see Appendix). 


To evaluate the impact of misalignment, we considered the 
particular setting in which the error due to registration is a 
pure shift. Although being a simplified case, this helps to 
understand the general algorithm performance as a function 
of the parameter p and the level of misalignment. In this 
particular case the translation error can be absorbed in a phase 
shift of the kernel. Although the weights will not change, since 
they only depend on the Fourier magnitude, the average in 0 
will be out-of-phase due to the misalignment of the images 
Vi. This will result in blur but also on image artifacts. 

We carried out a similar empirical analysis as before, where 
several thousands kernels were simulated and centered by forc¬ 
ing to have vanishing first moment. We introduce a Gaussian 
random 2d shift to each blurring kernel (i.e., the first moment 
of the kernel is shifted) with zero mean and standard deviation 
controlled by a parameter e to simulate the misalignment. 
Figure shows the algorithm performance when changing 
the amountof error in the registration (from 0-5 pixels in 
average) and the value of p controlling the FBA aggregaion. 
As the alignment error is more significant, the mean square 
error increases as expected. When the misalignment error is 
large, the best is to use low p values (close to arithmetic 
mean). On the other hand, when the misregistration error is 
low, large p values will produce the best performance in terms 
of reconstruction error (MSE). For shift errors in the order of 
1 pixel, the p value giving the minimum MSE is in p G [7, 20]. 
In general, the bias is always reduced with p while the variance 
is increased. This is a direct consequence of the fact that 
smoothness is reduced when increasing the p value (see in 
Fig. how the energy is concentrated in fewer images as 
p ^ oo, indicating less smoothing). 

D. Comparison to Classical Lucky Imaging Techniques 

Lucky imaging techniques, very popular when imaging 
through atmospheric turbulence, seek to select and average the 
sharpest images in a video. The Fourier weighting scheme can 
be seen as a generalization of the lucky imaging family. Lucky 
imaging selects (or weights more) frames/regions that are 
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(a) Examples of simulated kernels (exposure Tgxp) 



s = 0.01 s = 0.02 s = 0.04 s = 0.08 s = 0.16 
(b) Examples of noisy simulations (noise level s) 




(c) Noise level s 




(d) Number of images M 



Fig. 5. Bias-Variance tradeoff, (a) Simulated kernels due to hand tremor following (3- Each row shows a set of simulated kernels (left panel) for different 
exposures Texp = and the respective Fourier spectrum magnitude (bottom panel). The parameter Texp controls the amount of expected 

blur, (b) Simulated noise levels. Average algorithm performance with respect to p when changing (c) the amount of noise s in the input images, (d) the 
number of images in the burst M, and (e) the exposure time Texp. The rest of the parameters are set to M = 16, s = 0.04 and Texp = unless other 
specified. With short exposures, the arithmetic average (p = 0) produces the best mse since the images are not blurred. The bias does not depend on M, 
but the variance can be significantly reduced by taking more images (light accumulation procedure). Noise affects the bias and the variance terms (with the 
exception of p = 0 where the bias is unaffected). The mse plots show the existence of a minimum for p G [7, 30], indicating that the best is to balance a 
perfect average and a max pooling. 



(e) Registration error e 


Fig. 6. Impact of misalignment. We simulated shifts following a 2d Gaussian 
distribution of zero mean and standard deviation e. The average displacement 
is approximately 1.25e (shown in the table). The top-right plot shows the 
average algorithm performance (in MSE) with respect to p when changing 
the amount of registration error e. Misalignment deteriorates the algorithm 
performance. When the error is large (e > 2), low p values (aggregation is 
close to arithmetic average) produces the best mse, while when the registration 
error is low (e < 1), large p values produce lower reconstruction errors. Using 
a large p produces results that are less smoothed, so in general, variance is 
increased while bias is reduced. 


sharper but without paying attention to the characteristics of 
the blur. Thus, when dealing with camera shake, where frames 
are randomly blurred in different directions, classical lucky 
imaging techniques will have a suboptimal performance. In 
contrast, trying to detect the Fourier frequencies that were less 
affected by the blur and then build an image with them makes 
much more sense. This is what the Fourier Burst Accumulation 
seeks. 

Recently, Garrel et al. m introduced a selection scheme for 
astronomic images, based on the relative strength of signal 
for each frequency in the Fourier domain. Given a spatial 
frequency, only the Fourier values having largest magnitude 
are averaged. The user parameter is the percentage of largest 
frames to be averaged (typically ranging from 1%-10%). This 
method was developed for the particular case of of astro¬ 
nomical images, where astronomers capture videos having 
thousands of frames (for example 9000 frames in 171). Our 
algorithm is built on similar ideas, but we do not specify 
a constant number of frames to be averaged. We let the 
Fourier weighting scheme select the total contribution of 
each frame depending on the relative strength of the Fourier 
magnitudes. The authors showed that this generalized lucky 
imaging procedure produces astronomical images of higher 
resolution and better signal to noise ratio than traditional lucky 








































































Input frames crops 1-12 ordered from highest (left) to lowest (right) gradient energy 



Arithmetic average Sharpness-selectivity lfa 1 lfa 2 lfa 3 FBA p = 11 


Fig. 7. Comparison to lucky imaging techniques on a real data burst (building). The arithmetic average produces the best noise reduction but completely 
removes the details. The lucky imaging algorithm proposed in El] (sharpness-selectivity) slightly improves the arithmetic average. In this particular experiment 
there is blur in different directions and no single frame sharp in all directions. The sharpest frame detected (LFA 1) is still blurred and significantly noisier 
than the result given by the Fourier Burst Accumulation with p = 11 (final sharpening step has been disabled to do a fair comparison). As more sharp frames 
are averaged (lfa 2 and lfa 3) noise is reduced at the expense of blurring the final image. 


image fusion schemes, further stressing our findings. 

Traditional Lucky imaging techniques, are based on eval¬ 
uating the quality of a given frame. In astronomy, the most 
common sharpness measure is the intensity of the brightest 
spot, being a direct measure of the concentration of the 
system’s point spread function. However, this is not applicable 
in the context of classical photography. 

Haro et al. 1^ propose to locally estimate the level of 
sharpness using the integral of the energy of the image gradient 
in a surrounding region (i.e., the Dirichlet energy). If all the 
images in the series have similar noise level, the Dirichlet 
energy provides a direct way of ordering the images according 
to sharpness, (i.e., large Dirichlet energy implies sharpness). 

Let Vi i = ^ M, be a series of registered images of the 

same scene, the per-pixel Dirichlet energy weights are 

W^dirichlet(x) = [ |VVi(x)|2dx, (7) 

where is a block of (100 x 100) pixels around the pixel 
X. In practice, the Dirichlet weights vary poorly with camera 
shake blur, so although blurry images will contribute less to 
the final image, their contribution will be still significant. 

Joshi and Cohen ED propose to use a combination of 
sharpness and selectivity per-pixel weights to determine the 
contribution of each pixel to the restored image. The sharpness 
weight is built from the local intensity of the image Laplacian 
and it is pondered by a local selectivity term. The selectivity 
term enforces more noise reduction in smooth/fiat areas. The 
final sharpness-selectivity weights ar^ 

«^sWsel(^) = (8) 

where = max^]A«^(x)| sharpness measure 

and 7(x) = Al^^(^)l/max 3 , |af(x)| is the selectivity term con¬ 
trolled by a parameter A. We have denoted by v the average 
of all input frames. 

Tn the present analysis we did not consider the terms due to the resampling 
error nor the sensor dust that were originally included in the formulation 
presented in ED 


In all cases, the final image is computed as a per-pixel 
weighted average of the input images. 

Figure [7] shows a comparison of these two approaches 
to the proposed Fourier Burst Accumulation. We did not 
include a final sharpening step to faithfully compare all the 
approaches, as this last step could be included in all of them 
(see Section S- Since the weighted averaged image using 0 
did not show any difference with respect to the arithmetic 
average, we instead used the total Dirichlet energy to rank all 
the input images and then average only the top K (the sharpest 
ones). We named this method Lucky frame average (LFA) and 
tested different values of K = 1, 2,3. 

The weights given by ^ lead to an over-smoothed image 
with significant less noise. The sharpest frame detected (LFA, 
K = 1) is still blurred and noisier than the result given 
by the Fourier Burst Accumulation with p = 11. As more 
lucky frames are averaged (lfa 2 and LFA 3), more noise 
is eliminated at the expense of introducing blur in the final 
image. 



V. Algorithm Implementation 

The proposed burst restoration algorithm is built on three 
main blocks: Burst Registration, Fourier Burst Accumulation, 
and Noise Aware Sharpening as a post-processing. These are 
described in what follows. 

Burst Registration. There are several ways of registering 
images (see Ea for a survey). In this work, we use image cor¬ 
respondences to estimate the dominant homography relating 
every image of the burst and a reference image (the first one 
in the burst). The homography assumption is valid if the scene 
is planar (or far from the camera) or the viewpoint location is 
fixed, e.g., the camera only rotates around its optical center. 
Image correspondences are found using SIFT features 13611 and 
then filtered out through the ORSA algorithm l37]| . a variant 
of the so called RAN SAC method l3^ . To mitigate the effect 
of the camera shake blur we only detect SIFT features having 
a larger scale than cTmin = 1.8. 
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Recall that as in prior art, e.g., 1^ . the registration can be 
done with the gyroscope and accelerometer information from 
the camera. 


Fourier Burst Accumulation. Given the registered images 
we directly compute the corresponding Fourier trans¬ 
forms {vi}fLi. Since camera shake motion kernels have a 
small spatial support, their Fourier spectrum magnitudes vary 
very smoothly. Thus, \vi\ can be lowpass filtered before 
computing the weights, that is, \vi\ = Gcr|fii|, where is 
a Gaussian filter of standard deviation a. The strength of 
the low pass filter (controlled by the parameter a) should 
depend on the assumed motion kernel size (the smaller the 
kernel the more regular its Fourier spectrum magnitude). 
In our implementation we set a = where 

ks = 50 pixels and the image size is mhXrriw pixels. Although 
this low pass filter is important, the results are not too sensitive 
to the exact value of a (see Figure [^. 

The final Fourier burst aggregation is (note that the smooth¬ 
ing is only applied to the weights calculation) 


Urrj - 


f M 

E 

Vi=l 


Wi • Vi 


Wi = 




(9) 


The extension to color images is straightforward. The 
accumulation is done channel by channel using the same 
Fourier weights for all channels. The weights are computed 
by arithmetically averaging the Fourier magnitude of the 
channels before the low pass filtering. 


V ^ 
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Fig. 8. Impact of smoothing the Fourier weights. To eliminate image artifacts 
and noise, \vi \ are smoothed before computing the weights wi. Top row shows 
the results of the fba average (for the burst shown in Figure [7}, middle 
row the Fourier weights, and the bottom row a crop of the Fourier weights, 
when considering different levels of Gaussian smoothing (no smoothing, cr, 
3(7, 6a). The strength of the low pass filter is controlled by the parameter 
cr = min{mh,muj)/ks, where ks = 50 pixels and the image size is rrih x niw 
pixels. As shown in the left column (no smoothing), this filtering step is very 
important. It provides stabilization to the Fourier weights and also helps to 
remove noise. The results are stable for a large range of smoothing levels. 


are still available. The sharpening must contemplate that the 
reconstructed image may have some remaining noise. Thus, 
we first apply a denoising algorithm (we used the NLBAYES 
algorithm 1390. then on the filtered image we apply a 
Gaussian sharpening. To avoid removing fine details we finally 
add back a percentage of what has been removed during the 
denoising step. 

The complete method is detailed in Algorithm 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


Algorithm 1: Aggregation of Blurred Images 

Input: A series of images fii, i} 2 , • • •, of size 
m X n X c* . An integer value p. 


Output: The aggregated image Up 
w = zeros(m, n)\ Up = zeros(m, n, c); 
for image i in {1,..., n} do 
Burst Registration 

Mi = SIFT(i} 2 : set of corresponding points. 

Hi = ORSA(7lTy ; Hi dominant homography in Mi. 

Vi = Vi o Hi ; Image resampling. 


Fourier Burst Accumulation 
Vi = FFT(vi); 

~ c E^=l I’ Mean over color channels 

Wi = GtjWi ; Gaussian smoothing 

Up = Up • Vi ; Weighted Fourier accumulation 

W = W ^ W^\ 

Up = lFFT{Up/w); 

Noise Aware Sharpening (Optional) 

Up = DENOISE(lXp); 

= ‘lUp — GpUp\ Gaussian sharpening, p £ [1,3] 

Up = Up + 6(^Up — Up)’, Add a fraction of removed noise, (5 = 0.4 


HI c is the number of color channels, typically 3. The color channels 
of V are denoted by vf for j = 1,.. ., c. All the regular operations 
(e.g, -h, /, •) are point-wise. 


Memory and Complexity Analysis. Once the images are 
registered, the algorithm runs in 0{M • m • logm), where 
m = rrih x is the number of image pixels and M 
the number of images in the burst. The heaviest part of the 
algorithm is the computation of M FFTs, very suitable and 
popular in VLSI implementations. This is the reason why 
the method has a very low complexity. Regarding memory 
consumption, the algorithm does not need to access all the 
images simultaneously and can proceed in an online fashion. 
This keeps the memory requirements to only three buffers: one 
for the current image, one for the current average, and one for 
the current weights sum. 

VI. Experimental Results 


Noise Aware Sharpening. While the results of the Fourier 
burst accumulation are already very good, considering that 
the process so far has been computationally non-intensive, 
one can optionally apply a final sharpening step if resources 


We captured several handheld bursts with different number 
of images using a Canon 400D DSLR camera and the back 

^A variant of this is already available on camera phones, so we stay at the 
level of potential on-board implementations. 
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woods 13 imgs parking night 10 imgs booksheif 10 imgs 

ISO 1600, 1/8” ISO 1600, 1/3” ISO 100, 1/6” 

Canon 400D Canon 400D Canon 400D 



portrait 10 imgs 
ISO 800, 1/8” 
Canon 400D 


auvers 12 imgs 
ISO 400, 1/2”, iPad 


anthropologie [2^ 8 imgs 
ISO 100, 353 ms 


tequila [2^ 8 imgs 
ISO 100, 177 ms 


Fig. 9. Restoration of image bursts captured with two different cameras. The 
number of frames, the ISO sensitivity and the exposure time is indicated for 
each burst. Note that in the case of the iPad tablet, the exposure time is the one 
indicated by the app., and there is no guarantee that this is the real exposure 
time. Full images are available at the project’s website. 


The proposed method clearly outperforms 1^ in all the 
sequences. This algorithm introduces strong artifacts that de¬ 
graded the performance in most of the tested bursts. Tuning 
the parameters was not trivial since this algorithm relies on 
4 parameters that the authors have linked to a single one 
(named 7). We swept the parameter 7 to get the best possible 
performance. 

Our approach is conceptually similar to a regular align 
and average algorithm, but it produces significantly sharper 
images while keeping the noise reduction power of the average 
principle. In some cases with numerous images in the burst 
(e.g., see the parking night sequence), there might already 
be a relatively sharp image in the burst (lucky imaging). Our 
algorithm does not need to explicitly detect such “best” frame, 
and naturally uses the others to denoise the frequencies not 
containing image information but noise. 


camera of an iPad tablet. The full restored images and the 
details of the camera parameters are shown in Figure 
and Figure The photographs contain complex structure, 
texture, noise and saturated pixels, and were acquired under 
different lighting conditions. All the results were computed 
using p = 11. The full high resolution images are available at 
the project’s website 


A. Comparison to Multi-image Blind Deblurring 

Since this problem is typically addressed by multi-image 
blind deconvolution techniques, we selected two state-of-the- 
art algorithms for comparison m, 1^ . Both algorithms are 
built on variational formulations and estimate first the blurring 
kernels using all the frames in the burst and then do a step 
of multi-image non-blind deconvolution, requiring significant 
memory for normal operation. We used the code provided by 
the authors]^ The algorithms rely on parameters that were 
manually tuned to get the best possible results. We also 
compare to the simple align and average algorithm (which 
indeed is the particular case p = 0). 

Figures 10 and 12 show some crops of the restored 
images by all the algorithms. In addition, we show two 
input images for each burst: the best one in the burst and 
a typical one in the series. The proposed algorithm obtains 
similar or better results than the one by Zhang et al O, at 
a significantly lower computational and memory cost. Since 
this algorithm explicitly seeks to deconvolve the sequence, 
if the convolution model is not perfectly valid or there is 
misalignment, the restored image will have deconvolution 
artifacts. This is clearly observed in the bookshelf sequence 
where O produces a slightly sharper restoration but having 
ringing artifacts (see Jonquet book). Also, it is hard to read the 
word “Women” in the spine of the red book. Due to the strong 
assumed priors, O generally leads to very sharp images but 
it may also produce overshooting/ringing in some regions like 
in the brick wall (parking night). 


^ http://dev.ipol.im/~mdelbra/fba/ 

^ http://zoi.utia.cas.cz/files/fastMBD.zip 
https://drive.google.eom/file/d/0BzoBvkfRHe5bUF9jQlZsWXRYSkk/ 


B. Execution Time 

Once the images are registered, the proposed approach runs 
in only a few seconds in our Matlab experimental code, while 
O needs several hours for bursts of 8-10 images. Even if 
the estimation of the blurring kernels is done in a cropped 
version (i.e., 200 x 200 pixels region), the multi-image non¬ 
blind deconvolution step is very slow, taking several hours for 
6-8 megapixel images. 

C. Multi-image Non-blind Deconvolution 

Figure^] shows the algorithm results in two sequences pro¬ 
vided in 1^ . The algorithm proposed in 1^ uses gyroscope 
information present in new camera devices to register the burst 
and also to have an estimation of the local blurring kernels. 
Then a more expensive multi-image non-blind deconvolution 
algorithm is applied to recover the latent sharp image. Our 
algorithm produces similar results without explicitly solving 
any inverse problem nor using any information about the 
motion kernels. 

D. HDR imaging: Multi-Exposure Eusion 

In many situations, the dynamic range of the photographed 
scene is larger than the one of the camera’s image sensor. A 
popular solution to this problem is to capture several images 
with varying exposure settings and combine them to produce a 
single high dynamic range high-quality image (see e.g., iQi, 

ED). 

However, in dim light conditions, large exposure times are 
needed, leading to the presence of image blur when the images 
are captured without a tripod. This presents an additional 
challenge. A direct extension of the present work is to capture 
several bursts, each one covering a different exposure level. 
Then, each of the bursts is processed with the FBA procedure 
leading to a clean sharper representation of each burst. The 
obtained sharp images can then be merged to produce a high 
quality image using any existent exposure fusion algorithm. 

Figure shows the results of taking two image bursts with 
two different exposure times, and separately aggregating them 
using the proposed algorithm. We then applied the exposure 
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Typical Shot 


Best Shot 


Ahgn and average Sroubek & Milanfar (2^ Zhang et al. Proposed Method Proposed method 

(no final sharp.) 


Fig. 12. Real data burst deblurring results and comparison to state-of-the-art multi-image blind deconvolution algorithms (woods rows 1-3, parking night 
rows 4-6, bookshelf rows 7-8, sequences). 
























































Typical Shot Best Shot Align and average 



Sroubek & Milanfar Zhang et al. O Proposed method 

Fig. 10. Real data burst deblurring results and comparison with multi-image 
blind deconvolution methods (auvers). 


fusion algorithm of ll4^ to get a clean tone mapped image 
from the two burst representations. The fusion is much sharper 
and cleaner when using the Fourier Burst Accumulation for 
combining each burst than the one given by the arithmetic 
average or the best frames in each burst. 

VII. Conclusions 

We presented an algorithm to remove the camera shake blur 
in an image burst. The algorithm is built on the idea that each 
image in the burst is generally differently blurred; this being 
a consequence of the random nature of hand tremor. By doing 
a weighted average in the Fourier domain, we reconstruct 
an image combining the least attenuated frequencies in each 
frame. Experimental results showed that the reconstructed 
image is sharper and less noisy than the original ones. 

This algorithm has several advantages. First, it does not 
introduce typical ringing or overshooting artifacts present 
in most deconvolution algorithms. This is avoided by not 
formulating the deblurring problem as an inverse problem of 
deconvolution. The algorithm produces similar or better results 
than the state-of-the-art multi-image deconvolution while be¬ 
ing significantly faster and with lower memory footprint. 

We also presented a direct application of the Fourier Burst 
Accumulation algorithm to HDR imaging with a hand-held 


Typical Shot Best Shot Align and average 



Fig. 11. Real data burst deblurring results and comparison with multi-image 
blind deconvolution methods (portrait). 

camera. As a future work, we would like to incorporate a 
gyroscope registration technique, e.g., 1^ . to create a real¬ 
time system for removing camera shake in image bursts. 

A very related problem is how to determine the best capture 
Strategy. Giving a total exposure time, would it be more 
convenient to take several pictures with a short exposure (i.e., 
noisy) or only a few with a larger exposure time (i.e., blurred)? 
Variants of these questions have been previously tackled 1^ . 
1341 . l43l in the context of denoising / deconvolution tradeoff. 
We would like to explore this analysis using the Fourier Burst 
Accumulation principle. 

Appendix 

We consider that a burst is correctly aligned if each Vi 
satisfies Vi = u ki ^ rii, with the blurring kernel ki 
having vanishing first moment. That is, f /c^(x)x(ix = 0. This 
constraint on the kernel implies that the kernel does not drift 
the image u, so each Vi is aligned to u. 

An intuitive way to motivate this requirement is by analyz¬ 
ing the result of iteratively applying the blurring kernel a large 
number of times. Let /c(x) be a non-negative blurring kernel, 
with /j, = f k(x)xdx and = f k(x)(x — fi){x — fiYdx. 
Then, 




G {n/j., nTd) 
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building, 2x12 imgs, l/s” and 1/40”, iso 800, Canon 400D 



FBA - long exposure (LE) FBA - short exposure (SE) EBA - exposure fusion 




Fig. 14. HDR burst exposure fusion. Two different examples of fusion of two 12-image bursts captured with two different exposure levels. On top, the FBA 
average of each burst (left and middle) and the HDR fusion of these two images using da (right). Below, image crops of the best shot (sharpest, manually 
selected) in each burst and another typical in the series; the exposure fusion using the regular align and average to combine all the images of a burst, the 
exposure fusion using the best shot in each series, and the fusion using the proposed Fourier weighted average. Some of the crops have been rescaled to 
improve their contrast. 
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Input 1 Park & Levoy [2^ Proposed method 

Fig. 13. Restoration results with the data provided in (2^ (anthropologie 
and tequila sequences). 


where G (n/Lt, nS) is a Gaussian function with mean nii and 
variance nH. This is a direct consequence of the Central Limit 
Theorem. This Gaussian function can be decomposed into two 
different components: a centered Gaussian kernel (the low pass 
filter component) and a shifting kernel given by a Dirac delta 
function, that is, 

G {nil, nE) = G (0, nE) ^ Sn^j,- 

This means that iteratively applying n times the kernel k is 
(asymptotically) equivalent to applying a Gaussian blur with 
variance nE and then shifting the image an amount nfi. Thus, 
we can say that the original kernel k{x.) drifts the image 
“in average” an amount fi. Therefore, if we do not want the 
image to be shifted, the blurring kernel should have zero first 
moment. Following this argument all the simulated kernels 
were centered by forcing fi = f /c(x)x(ix = 0. 
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